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TRANSIENT  COMPRESSIBLE  FLOW  IN  A  PIPING  NETWORK: 

A  SOLUTION  METHOD  AND  COMPUTER  SIMULATION 

INTRODUCTION 

This  report  describes  work  conducted  under  the  NU SC  Independent 
Research  and  Independent  Exploratory  Development  project  titled 
"Transient  Compressible  Pipe  Flow  Simulation."  The  main  objectives  of 
this  project  were  to  develop  numerical  techniques  for  calculating 
transient  compressible  flow  of  a  fluid  through  complex  piping  networks 
and  to  create  a  user-oriented  computer  simulation  sufficiently  general  to 
be  applicable  to  a  wide  variety  of  networks.  The  end  result  of  this 
effort  is  a  computer  program  called  COMP,  which  has  been  implemented  at 
the  Naval  Underwater  Systems  Center  on  a  Digital  Equipment  Corporation 
VAX  11/780  computer. 

Several  simulations  for  adiabatic  flow  of  compressible  gas  in 
relatively  complex  piping  networks  have  been  developed  over  the  past  20 
years  (see  references  1  through  5).  Each,  however,  was  written  for  a 
specific  application  and  consequently  has  serious  drawbacks  when  applied 
to  other  networks.  Specifically,  the  simulations  are  not  user-oriented 
(making  it  difficult  to  define  and  set  up  the  network),  are  not  applicable 
to  choked  flow  across  throttling  devices  such  as  nozzles  and  valves,  and 
are  not  readily  applicable  to  transient  flow  where  a  combination  of 
conditions  (e.g.,  constant  pressure,  constant  volume,  or  variable  volume) 
may  exist  at  supply  and  receiving  tanks. 

The  simulation  described  in  this  report  overcomes  these  drawbacks. 

In  addition,  it  provides  a  new  technique  for  rapid  convergence  of  mass 
flow  rate  distributions  at  junctions  common  to  two  or  more  pipe  lines, 
can  account  for  any  ideal  gas,  and  employs  numerical  techniques  for  all 
calculations  requiring  iteration  or  integration.  It  can  be  expanded  to 
accommodate  practically  any  piping  network. 
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The  simulation  assumes  adiabatic  flow  in  each  line.  It  accounts  for 
irreversibilities  encountered  in  throttling  devices  and  for  frictional 
effects  in  the  pipe  and  fittings  by  treating  them  as  equivalent  pipe 
lengths.  Figure  1  shows  the  network  used  in  the  simulation.  Circled 
numbers  in  the  figure  indicate  tanks,  uncircled  numbers  indicate  pipe 
sections,  and  circled  letters  indicate  Junctions.  Pipe  sections  connect 
two  junctions  or  a  tank  and  a  junction.  Each  pipe  section  can  consist  of 
up  to  10  pipe  lengths  arranged  in  series.  Throttling  devices  can  be 
substituted  for  pipe  lengths,  but  cannot  be  the  first  or  last  lengths  in 
a  section. 
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DISCUSSION 


THEORY 

The  general  approach  used  here  in  solving  the  transient  compressible 
flow  problem  in  a  piping  network  is  similar  to  that  used  in  reference  1. 
Namely,  adiabatic  frictional  flow  is  assumed  to  exist  in  each  line  of  the 
network  and  the  system's  mass  flow  rate  is  determined  by  balancing  mass 
flow  rate  and  stagnation  pressure  at  a  single  point  (called  the  common 
point),  which  is  subject  to  the  total  system  mass  flow  rate.  For  the 
initial  time  increment  in  this  transient  flow  problem,  calculations  first 
proceed  in  a  forward  direction  from  the  most  remote  supply  tank  to  the 
common  point,  and  then  in  a  backward  direction  from  the  most  remote 
receiving  tank  to  the  common  point.  Stagnation  pressure  values  at  the 
common  point  as  calculated  by  these  forward  and  backward  computations  or 
"passes"  are  then  compared.  If  they  are  not  approximately  equal,  the 
mass  flow  rate  calculations  undergo  repeated  iterations  until  they 
balance.  Once  convergence  of  stagnation  pressures  at  the  common  point 
occurs,  the  simulation  employs  a  Runge-Kutta  integration  scheme  to  set  up 
tank  conditions  for  the  beginning  of  the  next  time  increment,  whereupon 
the  calculation  process  for  stagnation  pressures  at  the  common  point  is 
repeated . 

The  equations  used  to  calculate  compressible  adiabatic  flow  are  taken 
from  reference  6  and  are  not  repeated  here.  The  convergence  scheme, 
however,  was  developed  along  with  the  simulation  itself.  The  scheme  is 
programmed  as  one  of  the  simulation's  main  subroutines  (subroutine  RATIO) 
and  is  discussed  along  with  other  subroutines  later. 


ASSUMPTIONS 


Assumptions  made  in  conducting  this  analysis  are  as  follows: 

1.  Flow  is  one-dimensional,  quasi-steady,  and  adiabatic. 

2.  Ratio  of  specific  heats  is  constant. 

3-  Ideal  fluid  compressibility  factor,  Z,  is  constant  and  equal  to 
unity. 

4.  Friction  factor,  f,  for  any  length  of  pipe  i3  equal  to  the 
average  of  the  inlet  and  outlet  steady-state  friction  factors. 

5.  Minor  losses  due  to  elbows,  reducers,  etc.,  are  introduced  as 
equivalent  pipe  lengths. 

6.  No  reverse  flow  exists  in  any  line. 

7.  Isentropic  stagnation  pressures  are  balanced  at  each  junction  of 
two  or  more  pipe  lines. 

8.  Isothermal  thermodynamic  processes  exist  in  each  tank. 


STRATEGY 

The  junctions  between  pipe  sections  are  assumed  to  be  isentropic, 
since  any  frictional  effects  due  to  the  junctions  may  be  compensated  for 
by  adding  equivalent  pipe  lengths  to  the  network.  The  isentropic 
condition  at  junctions  forces  the  conclusion  that  there  must  be  a  single 
common  stagnation  pressure  for  all  pipe  lengths  common  to  a  junction. 
This  assumption  is  critical  to  the  analysis  because  it  leads  to  the 
condition  that  the  stagnation  pressure  is  the  property  to  be  iterated 
upon  at  each  junction  in  order  to  determine  mass  flow  rate  distribution. 


One  important  exception  to  this  condition  is  evident  when  throttling 
or  choking  occurs  in  a  pipe  terminating  at  a  junction.  The  stagnation 
pressure  calculated  at  the  end  of  this  pipe  may  be  greater  than  the 
junction  stagnation  pressure,  with  the  pressure  difference  being  lost  in 
the  choked  pipe  as  the  flow  enters  the  junction. 

The  overall  approach  to  calculations  in  the  forward  direction  is  to 

assume  a  mass  flow  rate  in  the  most  remote  and  always  present  supply  tank. 

This  is  tank  1  and  mass  flow  rate  in  figure  1.  Throughout  the 

transient  there  must  be  some  absolute  flow  from  this  tank  since  it 

controls  the  overall  calculation  process.  Calculations  proceed  to  the 

first  junction  (junction  A  in  figure  1),  where  a  stagnation  pressure  at 

the  end  of  pipe  section  1  is  calculated.  The  program  next  assumes  a 

mass  flow  rate  from  tank  2  into  pipe  section  2  (m2)  ,  and  calculates  a 

stagnation  pressure  at  junction  A  from  section  2.  If  this  stagnation 

pressure  is  within  1  psi  of  that  calculated  from  section  1,  calculations 

proceed  to  section  4;  if  it  is  not,  adjustment  is  made  to  m2>  via  the 

convergence  scheme  (subroutine  RATIO),  while  m^  ^3  held  at  the  initial 

value.  As  soon  as  agreement  is  reached  (usually  2  to  5  iterations  on 
• 

m2),  an  arithmetic  average  of  stagnation  pressures  at  junction  A  is 
made  and  calculations  continue. 

•  •  • 

Mass  flow  rate  then  taken  to  be  the  sum  of  and  m2. 

The  starting  pressure  for  section  4  is  taken  as  the  junction  stagnation 
pressure  obtained  from  sections  1  and  2.  A  stagnation  pressure  at 
junction  B  from  section  4  is  then  calculated.  The  iterative  procedure 
used  to  calculate  the  tank  2  flow  is  next  employed  again  for  tank  3  to 
obtain  a  value  for  m^  and  a  stagnation  pressure  at  junction  B.  This 
completes  the  forward  pass,  since  junction  B  is  the  common  point  where 
the  forward  and  backward  calculations  meet. 
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The  backward  pass  begins  at  the  receiving  tanks  and  proceeds  toward 
the  common  point,  junction  B.  It  is  assumed  that  whatever  mass  flowed 
into  the  system  from  the  supply  tanks  will  flow  out  of  the  system  into 
the  receiving  tanks  during  any  particular  time  step.  To  ensure  that  this 
condition  is  met,  the  program  takes  a  percentage  of  the  mass  flowing  into 
junction  B  from  the  forward  calculations,  and  uses  that  as  the  mass  flow 
in  section  9*  The  remainder  goes  to  section  10.  The  flow  in  section  9 
is  then  divided  between  sections  11  and  12  (initially  50  percent  in 
each).  Calculations  proceed  to  determine  the  stagnation  pressure  at 
junction  F  that  would  be  required  to  produce  that  estimate  of  m.,  . 
Stagnation  pressure  at  junction  F  from  section  12  is  calculated  next  and 
compared  to  the  pressure  from  section  11.  If  the  two  do  not  agree,  a 
mass  flow  rate  adjustment  is  again  made  via  subroutine  RATIO. 

Once  agreement  is  reached,  m^  j_s  used  to  back-calculate  the 
stagnation  pressure  required  at  junction  E  for  a  flow  rate  of  m^.  ^he 
program  proceeds  identically  in  sections  13  and  14  to  junction  E.  If  the 
stagnation  pressure  at  junction  E  from  section  9  does  not  match  the  one 

from  section  10  within  limits  set  into  the  program  (presently  1  psi),  then 

•  • 

subroutine  RATIO  adjusts  mg  and  mj_Q  until  agreement  is  reached.  The 
•  • 

correct  mg  and  m-^Q  values  are  rapidly  obtained,  resulting  in  a 

calculated  stagnation  pressure  at  junction  E  from  the  backward 

•  •  • 
calculation.  Note  that  if  mg  and  m^  are  adjusted,  then  m-j^ 

through  must  also  be  adjusted  and  all  backward  pass  calculations 
repeated  from  the  tanks. 

Calculations  proceed  to  the  beginning  of  section  8  where  stagnation 
pressure  at  junction  D  is  obtained.  The  flow  then  splits  into  the 
parallel  branches  of  sections  6  and  7.  Under  the  initial  assumption  that 
the  flow  is  equally  split  between  sections  6  and  7,  the  pressure  at 
junction  C  based  on  each  section's  contribution  is  then  determined. 
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If  agreement  is  achieved  within  the  set  limits  of  1  psi,  the  parallel 
branch  calculations  are  complete  and  the  simulation  calculates  the 
conditions  for  section  5.  If  there  is  no  agreement,  subroutine  RATIO 

adjusts  the  percentage  of  flow  going  to  section  6  (which  then  gives  a  new 

•  •  • 

flow  rate  in  section  7  since  m^  -  ^  _  mg)  until  convergence  is 
achieved.  Stagnation  pressure  at  junction  B  is  then  calculated.  This 
completes  the  backward  pass,  since  junction  B  is  the  common  point. 


The  overall  convergence  scheme  of  subroutine  RATIO  is  once  again 
accessed.  It  now  compares  the  stagnation  pressure  at  junction  B  from  the 
forward  pass  with  stagnation  pressure  from  the  backward  pass.  If  they 
match  (within  1  psi)  the  network  is  solved  for  that  time  step.  If  they 
do  not  match,  an  adjustment  to  m^  j.3  made  via  subroutine  RATIO  and  the 

forward  and  backward  pass  calculations  are  again  repeated  with  the  new 

•  • 

value  for  m^ .  Iterations  on  m^  continue  until  stagnation  pressure 
values  at  common  point  B  as  calculated  from  the  forward  and  backward 
passes  agree.  At  that  time,  all  distributed  mass  flow  rates  and  junction 
pressure  values  are  correct,  values  for  that  particular  time  increment 
are  output  to  the  output  file,  and  calculations  for  subsequent  time 
increments  begirt.  Initial  conditions  used  for  the  subsequent  time 
increment  are  determined  by  a  Runge-Kutta  integration. 


RUNGE-KUTTA  INTEGRATION  SCHEME 

The  accuracy  of  the  program  in  predicting  the  transient  behavior  of 
the  network  depends,  in  part,  upon  the  technique  selected  to  update  the 
conditions  within  the  supply  and  receiving  tanks.  The  experience  of  the 
authors  with  different  numerical  integration  schemes  led  to  the  selection 
of  a  fourth-order,  Runge-Kutta  technique,  which  is  ideally  suited  to  this 
particular  application. 
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The  technique  is  straightforward.  The  initial  condition  of  each  tank 
with  respect  to  volume,  temperature,  pressure,  and  type  of  gas  being  used 
is  known.  The  initial  mass  of  gas  contained  within  each  tank  is  first 
calculated  by  the  perfect  gas  law,  the  forward  and  backward  passes  are 
run,  and  the  overall  convergence  criteria  are  applied  at  the  common 
point.  Once  convergence  at  the  common  point  is  achieved,  the  mass  flow 
rate  from  (or  into)  each  tank,  as  well  as  that  through  all  the  pipes,  and 
across  each  junction  is  known.  With  the  initial  mass  within  a  tank  and 
mass  flow  rate  known,  it  is  possible  to  calculate  a  new  mass  within  the 
tank  by  means  of  a  Runge-Kutta  integration.  Figure  2  illustrates  this 
concept.  Once  the  new  mass  is  found,  the  perfect  gas  law  is  once  again 
used  to  obtain  new  pressures  under  the  assumption  of  isothermal  tank 
processes.  Having  new  pressures  at  the  beginning  of  the  next  time 
increment  allows  the  calculation  process  to  begin  again. 
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Figure  2.  Runge-Kutta  Integration  Scheme 


In  the  event  that  the  tank  has  a  variable  volume,  the  user  may  employ 
a  differential  equation  to  describe  the  volume  rate  of  change  of  the  tank 
with  respect  to  time.  The  differential  equation  can  be  non-linear,  but 
must  be  ordinary.  In  running  the  simulation,  the  user  inserts  the 
equations  into  the  FORTRAN  main  program  in  a  section  clearly  labeled  with 
comment  statements.  An  example  would  be  the  following  equation: 

dV2/dt  =  2  x  Time  x  Time  +32.0 

which  would  be  inserted  into  the  program  as 

F(2)  =  2  *  TIME  *  TIME  +  32-0 

in  the  section  marked  for  tank  2.  It  would  also  be  necessary  to 
initialize  the  volume  of  tank  2  in  this  example  at  time  =  0.  The  FORTRAN 
statement  is  Y(2)  =  32.0,  which  would  indicate  an  initial  volume  of  32 
cubic  feet. 

CONVERGENCE  TECHNIQUE 

The  assumptions  that  the  flow  is  adiabatic  and  frictional  provide  no 
explicit  equation  for  pressure  drop  in  terms  of  a  mass  flow  rate  through 
a  pipe,  as  would  be  the  case  for  isothermal  pipe  flow.  Instead,  the 
problem  involves  a  set  of  very  non-linear  algebraic  equations  that  must 
be  solved  simultaneously.  The  problem  experienced  when  the  wrong  mass 
flow  rate  is  initially  selected  is  that  the  stagnation  pressures  do  not 
match  at  the  junctions.  For  example,  consider  section  2,  with  m^ 
flowing  to  junction  A.  At  junction  A  a  value  for  stagnation  pressure 
from  section  1  has  been  obtained  and  it  is  now  necessary  to  arrive  at  the 
correct  value  of  m^  in  order  to  have  the  junction  A  stagnation  pressure 
from  section  2  match  the  pressure  from  section  1.  In  short  m^  must  be 
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made  to  converge  to  the  correct  value.  Since  the  equations  being  used  to 
calculate  junction  pressures  are  non-linear,  a  technique  had  to  be 
devised  to  cause  m 2  to  converge  rapidly.  Indeed,  one  of  the  main 
contributions  of  this  study  is  the  method  that  was  developed  to  cause  the 
mass  flow  rates  to  converge.  The  method  is  extremely  fast  and  general. 
For  example,  with  an  initial  estimate  of  a 2  0f  q.1  Ibm/s  when  the 
actual  converged  value  should  be  100.0  lbm/s,  the  method  used  here  will 
yield  the  correct  value  in  only  six  iterations. 


The  method  is  as  follows:  an  initial  value  of  m2  i3  estimated  from 
which  a  stagnation  pressure  at  junction  A  is  calculated.  If  the  estimate 
does  not  result  in  convergence  it  must  be  adjusted  and  the  procedure 
repeated.  If  the  stagnation  pressures  from  sections  1  and  2  are  far 
apart  at  junction  A,  then  a  relatively  large  adjustment  must  be  made,  and 
the  direction  in  which  the  adjustment  must  be  made  is  known.  If  the 
value  calculated  for  section  2  was  too  low  at  the  junction,  then  m2 
must  be  reduced  so  there  is  less  pressure  drop  through  section  2,  thus 
increasing  the  section  2  stagnation  pressure  value  at  the  junction.  If 
the  calculation  for  section  2  was  too  high  at  the  junction,  then  m2 
must  be  increased  in  order  to  lower  the  section  2  pressure  value.  It  is 
not  desirable  to  simply  add  or  subtract  a  fixed  amount  to  m2  because  it 
is  not  possible  to  know  ahead  of  time  how  large  m2  should  be. 

For  efficient  convergence,  a  method  was  developed  that  is  sensitive 
not  only  to  the  size  (and  sign)  of  the  pressure  difference  at  junction 
A,  but  to  the  size  of  m2  as  well.  The  method  uses  the  exponential 
function : 


ynew  =  ^old^  exP  (Constant  *  x) 

where  y  =  m2  and  x  =  junction  A  stagnation  pressure  from  section  2 
minus  the  stagnation  pressure  from  section  1. 
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The  choice  of  constant  will  depend  on  how  quickly  convergence  is 
desired  and  how  sensitive  x  is  to  small  changes  in  y.  The  constant  value 
in* the  simulation  as  it  currently  exists  is  good  for  most  situations,  but 
may  have  to  be  changed  for  very  small  pipe  diameters.  Tables  la  and  lb 
illustrate  convergence  for  two  very  different  size  y-values  when  the 
constant  equals  0.00005* 

Table  la.  Small  y-Value  Convergence 


+  1000.0 
+  50.0 
-  10.0 
+  1.5 
-  0.2 


0.1 

0.1051271 
0.1053903 
0.1053376 
0. 1053*155 


Convergence 


0.1051271 

0.1053903 

0.1053376 

0.1053455 

0.1053444 


Table  lb.  Large  y-Value  Convergence 


+  1005.0 

85.0 

89.380385 

+  500.0 

89.380385 

91.643061 

300.0 

91.643061 

93.02806 8 

+  60.0 

93.028068 

93-307571 

-  5.0 

93-307571 

93-284247 

+  0.2 

93.284247 

Convergence 

93.285180 

This  convergence  method  is  also  used  to  adjust  the  ratio  of  mass  flow 
rates  at  a  junction  having  two  pipe  sections  leaving  or  entering  from  it. 
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At  junction  C,  for  example,  RATI0(6)  can  be  defined  as  0 tal ' 

Then,  by  adjusting  RATIQ(6)  using  two  equations  (see  figure  3)  it  is 
possible  to  split  the  flow  by  the  appropriate  amount: 

a  =  RATI0C6) 

new 

=  RATI0(6)old  expCC^x)  (for  >  0  and  x  >  0) 

=  1  -  exp(-C2X)  +  RATI0(6)old  *  exp(-C2x)  (for  C2  >  0  and  x  <  0) 


RATIO  (6) 


Figure  3*  Ratios  of  Mass  Flow  Rates 

The  following  summary  shows  where  in  the  piping  network  the  convergence 
scheme  of  subroutine  RATIO  is  applied,  as  well  as  the  the  order  of  use: 


1.  At  junction  A  to  adjust  m2 

2.  At  junction  B  to  adjust  m, 

3.  At  junction  F  to  adjust 

U.  At  junction  G  to  adjust 

5*  At  junction  E  to  adjust 

6.  At  junction  C  to  adjust 

7.  At  junction  B  for  overall  convergence,  to  adjust  .'n1 

Uses  1  and  2,  above,  are  for  the  forward  pass;  uses  3  through  6  are  for 


the  backward  pass. 


THROTTLING  DEVICE;  FORWARD  PASS 


Throttling  devices  control,  monitor,  or  limit  flow  rate  in  a  pipe. 
Examples  are  valves,  orifices,  and  nozzles.  The  effectiveness  and/or 
discharge  coefficient  of  such  devices  may  be  a  constant  or  may  be  a 
function  of  time,  Reynolds  number,  or  some  other  parameter.  Throttling 
devices  are  represented  schematically  as  shown  in  figure  4. 


Figure  4,  Throttling  Device  in  Pipe  Section 

The  equations  used  in  the  program  to  define  the  mass  flow  rate 
through  a  throttling  device  are  taken  from  reference  7.  These  equations, 
which  were  derived  for  orifices  with  pipe  taps,  are  as  follows: 


P  -  ? 

Y  ■  1  -  [0.333  +  1.145  (fi2  +  0.7a5  «.  12a13)]  1  p- 
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m  =  C0YA[2gP1(P1  -  P2)]“ 


(unchoked  condition) 


3*  1 1/2  v  /  9  \2/(Y-l )]  1/2 

i  si.  1  7TT  (my  i  (cl,oke<i  condici 
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P  =  P  I - 

2  crit  1  crit  l  y 


P  ,  :  P  ,  +  P  -  P, 
o2  ol  2  1 
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where: 


A  =  throttling  device  reference  area  (ft2) 

CD  =  discharge  coefficient 
g  =  acceleration  of  gravity  (ft/s2) 
m  =  mass  flow  rate  (lbm/s) 

=  upstream  pressure  (lb/ft2) 

=  downstream  pressure  (lb/ft2) 

R  =  gas  constant  [ (ft-lbf ) /( lbm-°R) ]  (53*3  for  air) 

T  =  temperature  (°R) 

S  =  diameter  ratio  (throttling  device/upstream  pipe) 
p  =  weight  density  (lb/ft^) 
y  =  specific  heat  ratio 

and  the  subscript  "o"  specifies  stagnation  properties  and  "crit"  specifies 
properties  at  critical  (throttled  or  choked)  flow. 

The  conditions  upstream  of  the  throttling  device  are  known  prior  to  the 
point  in  the  computer  program  at  which  the  throttling  device  calculation 
takes  place.  These  known  conditions  include  values  for  m,  (  ancj 

T01.  The  calculation  proceeds  as  follows: 

1.  The  critical  mass  flow  rate  (®crj_t)  is  calculated  via  equation  (3)< 

•  • 

2.  If  m  is  greater  than  ®crit,  then  the  estimate  of  flow  rate  from 
tank  1  is  too  high  and  is  reduced  by  ®orit/m.  Calculations 
begin  again  from  tank  1. 
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j*  ii  m  =  mcr^t,  then  the  flow  is  choked  and  will  pass  m.  For  this 
case  downstream  properties  are  calculated  via  equations  (4)  and  (5) 

•  •  • 

4.  If  m  is  less  than  ®cr^t  the  flow  is  not  choked  and  will  pass  m. 

Downstream  properties  are  calculated  next  via  equations  (1),  (2), 

and  ( 3 ) . 


THROTTLING  DEVICE;  BACKWARD  PASS 

When  a  throttling  device  is  encountered  in  the  backward  pass 

calculations,  downstream  conditions  (point  2  in  figure  4)  will  be  known 

(i.e.,  m,  Pq2>  and  To2)  •  T*ie  ma-n  difference  between  the  forward 

and  backward  pass  throttling  device  calculations  is  that  the  mass  flow  rate 

•  • 

is  not  iterated  in  the  backward  pass  when  m  is  greater  than  mQ .  Instead, 
the  upstream  pressure  is  adjusted  to  accommodate  the  required  mass  flow 
rate.  The  specific  procedure  is  as  follows: 

1*  Downstream  static  pressure,  j.s  calculated. 

2.  Critical  upstream  static  pressure,  ifc,  iS  calculated  via 

equation  (4). 

3.  Unchoked  upstream  static  pressure,  ?1  unchoked ,  is  calculated  via 
equations  (1)  and  (2).  Note  that  this  is  the  upstream  pressure 
required  to  produce  a  mass  flow  rate  across  the  device  equal  to  m. 

U*  If  Pl  unchoked  is  less  than  or  equal  to  Px  crit,  then 

unchoked  i3  the  upstream  static  pressure,  Pj_ ,  and  the 
upstream  stagnation  pressure  is  calculated  via  equation  (5). 
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5*  If  pi  unchoked  is  Sweater  than  Px  crit,  then  the  flow  across 
the  device  in  choked  and  upstream  stagnation  pressure,  PQj_ ,  is 
calculated  via  equation  (4).  In  this  case,  when  total  convergence 
for  the  entire  network  is  satisfied  all  the  properties  downstream 
of  the  throttling  device  will  have  to  be  recalculated  in  a  forward 
pass  fashion.  The  reason  for  this  is  that  the  mass  flow  rate  is  so 
high  in  this  line  that  there  is  a  loss  of  stagnation  pressure  from 
the  end  of  the  line  into  the  tank. 


This  concludes  the  calculations  across  a  throttling  device.  Although 
most  of  the  calculations  are  straightforward,  calculating  the  unchoked 
upstream  static  pressure,  P^  unchoked ’  3teP  3  can  be  difficult.  Here, 
this  was  overcome  by  means  of  a  Newton-Raphson  technique  (reference  8). 
This  technique  was  implemented  by  obtaining  a  single  algebraic  equation  at 
point  1  for  p^  as  a  function  of  m,  effective  area,  discharge  coefficient, 
specific  heat  ratio,  and  the  static  pressure  at  point  2.  This  equation, 
which  was  solved  by  the  Newton-Raphson  technique  for  p^,  is  as  follows: 


m 


C*  [2g(p  2RT  -  G  -  p.P.)]1/2 
I  o  12 


=  1 


-  —  (i  -  — — — — ) 

Y  \  (p.RT  -  r  )  ' 


(plRTo  - 


where 


C*  is  the  product  of  discharge  coefficient  (C^)  and  the  effective 
area  of  the  throttling  device, 


> 


and 


D  = 


’2 

m 
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Next,  an  equation  was  derived  to  calculate  the  static  pressure  at 
point  1  in  terms  of  p . . 


P  =  P  RT  -  — 

1  1  o  2gY  \  Ax  /  P]_ 


which  is  the  value  for  unchoked  static  pressure  at  point  1. 


PSOGRAM  STRUCTURE 


MAIN  PROGRAM 


The  main  program  (program  MAIN)  initializes  all  constants  and  reads 
all  input  data  from  an  input  file.  Next,  the  initial  estimate  for  the 
controlled  flo”  rate  from  tank  1  is  specified  via  the  terminal.  Once  all 
the  data  are  fed  into  the  computer,  program  MAIN  calls  various 
subroutines  to  calculate  the  pressures  and  flows  at  the  pipe  junctions 
throughout  the  system. 


Program  MAIN  calls  subroutine  CTANK  to  calculate  pressures  in  all 
pipes  leading  away  from  the  input  tanks.  Subroutine  REVERSE  is  called 
upon  to  calculate  the  pressures  and  flow  rates  from  the  receiving  tanks 
back  to  junction  D.  Once  the  flow  is  known  at  junction  D,  subroutine 
REVERSE  calculates  all  necessary  pressures  and  flow  rates  in  the  parallel 
branch.  REVERSE  then  evaluates  the  pressure  and  flow  rate  in  the  last 
pipe  (section  5)  of  the  backward  pass  to  the  common  point,  junction  B. 

Program  MAIN  then  checks  the  difference  in  stagnation  pressure  at  the 
common  point.  If  the  difference  is  large,  MAIN  will  change  the  control 
flow  rate  via  subroutine  RATIO,  and  will  repeat  the  entire  calculation 


procedure  until  agreement  is  reached  and  state  conditions  at  all 
junctions  of  the  flow  network  are  satisfied.  MAIN  will  print  the 
results,  increment  the  time  step,  and  repeat  the  procedure. 

Appendix  A  presents  flowcharts  of  the  main  program  and  major 
subroutines.  A  brief  description  of  the  subroutines  follows. 

MAIN  SUBROUTINES 

Subroutine  PRES  calculates  the  end  pressure  for  a  pipe  length  in  the 
forward  pass  under  the  assumption  that  specific  heat  ratio,  gas  constant 
mass  flow  rate,  stagnation  temperature,  stagnation  pressure,  pipe 
dimension,  and  pipe  length  are  known. 

Subroutine  TANK PR  calculates  tank  stagnation  pressure  at  the  end  of 
each  time  step  using  a  Runge-Kutta  integration. 

Subroutine  BKVALVE  performs  throttling  device  and  pipe  flow 
calculations  for  the  backward  pass. 

Subroutine  FVALVE  performs  throttling  device  calculations  for  the 
forward  pass. 

Subroutine  CTANK  calculates  end  pressures  and  flow  rates  of  pipes 
leaving  the  tanks  in  the  forward  pass  network  and  in  pipe  section  4. 

Subroutine  RATIO,  the  convergence  subroutine,  adjusts  the  flow  rate, 
using  an  exponential  function,  until  the  correct  value  for  the  flow  rate 
is  found. 


Subroutine  REVERSE  controls  the  backward  direction  calculations, 
which  calculate  end  pressure  and  flow  rate  for  all  junctions  between  the 
common  point  and  the  receiving  tanks.  It  also  calculates  all  pipe  length 
pressures  and  flow  rates  in  the  backward  pass. 

REAL  FUNCTION  SUBROUTINES 

Real  function  FRIC  determines  the  Darcy  friction  factor  at  a  point  in 
the  pipe  where  the  Mach  number,  stagnation  temperature,  stagnation 
pressure,  specific  heat  ratio,  gas  constant,  mass  flow  rate,  and  pipe 
diameter  are  known. 

Real  function  DYNVIS  uses  Sutherland's  equation  (reference  9)  to 
calculate  the  dynamic  viscosity. 

Real  function  MALPHA  calculates  the  Mach  number  at  the  end  of  the 
pipe  length  (outlet  condition),  using  fL*/D  ratios: 

7lvd  ,  (a»/B>inlet  .  (7LVD)outl8t 

where  L*  is  the  length  of  duct  required  to  develop  a  flow  from  the  Mach 
number  at  the  position  under  consideration  to  the  sonic  point,  f  is  the 
average  friction  factor,  and  D  is  the  pipe  diameter. 


Real  function  BETA  solves  the  equation  that  relates  Mach  number  to 
friction  given  a  known  Mach  number,  M,  and  specific  heat  ratio,  Y: 


( Y  ♦  1)M 


fL*/D  =  in 

YM*  2  +  (Y  -  1)M 


a 


Real  function  FF  is  a  digital  equivalent  of  the  Moody  chart  for 
friction  factor  as  a  function  of  pipe  roughness  and  Reynolds  number. 
(Units  are  in  feet  with  roughness  =  0.00015  for  commercial  steel  pipe.) 

Real  function  MARAT  calculates  the  Mach  number  at  the  pipe  inlet. 

Real  function  TEMP  calculates  the  temperature  at  any  point  given  the 
Mach  number  and  specific  heat  ratio. 

Real  function  CVEL  calculates  the  speed  of  sound  for  a  specified 
temperature. 

Real  function  VEL  calculates  the  velocity  of  the  fluid  at  any  point 
in  the  pipe  given  the  Mach  number. 

Real  function  RHO  calculates  the  density  of  the  fluid  at  any  point  in 
the  pipe  given  the  Mach  number  and  velocity. 

Real  function  REND  calculates  the  Reynolds  number  at  any  point  given 
the  Mach  number,  temperature,  density,  and  velocity. 


DEFINING  THE  NETWORK 


ORGANIZATION  OF  INPUT  DATA 


Input  data  for  the  compressible  pipe  flow  computer  program  is  set  by 
the  user  via  a  computer  file.  Figure  1  is  the  model  for  the  input  data, 
and  networks  up  to  the  complexity  shown  in  figure  1  can  be  simulated. 
Tank  1  must  always  be  present,  but  any  of  the  other  tanks  can  be  deleted 
as  long  as  at  least  one  receiving  tank  is  available.  Pipe  sections  can 
also  be  deleted  as  long  as  continuity  in  the  flow  path  is  maintained. 
Each  pipe  section  can  be  defined  by  up  to  10  individual  pipe  lengths, 
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each  with  its  corresponding  diameter.  Throttling  devices  are  substituted 
for  pipe  lengths,  but  they  cannot  be  located  as  the  first  or  last  lengths 
of  a  pipe  section.  This  is  of  no  serious  consequence,  since  if  a 
throttling  device  located  at  the  end  of  a  pipe  section  in  a  real  network 
must  be  simulated,  a  relatively  short  pipe  length  of  large  diameter 
(which  would  have  a  negligible  effect  on  the  calculations )  can  be  input 
in  the  simulation  after  the  throttling  device.  The  organization  of  the 
input  file  (see  appendix  B  for  a  sample  input  data  file)  is  as  follows: 

Line  1  -  Identification  line  for  the  problem  as  defined  in  columns  1 
to  110  of  the  data  file. 

Line  2  -  Stagnation  temperature  (°R)  in  the  network.  The  FORTRAN 
format  is  F10.4. 

Line  3  -  Initial  stagnation  pressure  (lb/ft^)  in  each  tank.  The 
format  for  this  line  is  7F10.4  and  is  in  the  following  tank  order:  1,  2, 
3,  11,  12,  13,  14.  If  a  tank  is  not  present,  stagnation  pressure  of  0.0 
is  input . 

Line  4  -  Initial  volume  (ft^)  of  each  tank.  The  format  for  this 
line  is  7F10.4.  Tank  order  is  the  same  as  for  line  3-  If  a  tank  is  not 
present,  volume  0.0  is  input. 

Line  5  -  Flags  specifying  which  tanks  are  in  the  system.  Foiuiat  is 
714.  Tank  order  is  the  same  as  for  line  3*  An  integer  value  of  1  is 
used  if  a  tank  is  in  the  system;  2  if  a  tank  is  not  in  the  system. 

Line  6  -  Flags  to  specify  the  conditions  in  each  tank.  The  format 
for  this  line  is  714  with  the  tank  order  the  same  as  for  line  3.  An 
integer  value  of  0  is  used  if  a  tank  is  not  present;  1  for  constant 
pressure;  2  for  constant  volume;  3  for  variable  volume  defined  by 
user-supplied  differential  equation. 
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Lines  7-8  -  The  number  of  different  pipe  lengths  in  each  pipe 
section.  Maximum  is  10.  Format  is  714.  Information  is  given  in  the 
following  order: 

Line  7:  Pipe  sections  no.  1,  2,  3.  4,  5,  6,  7 

Line  8:  Pipe  sections  no.  8,  9,  10,  11,  12,  13,  14 

If  a  pipe  section  is  not  present  in  the  network  an  integer  value  of  0  is 
input.  Each  throttling  device  counts  as  a  single  pipe  length. 

Lines  9-*-  -  Pipe  length  and  diameter  data  (ft)  for  each  pipe  listed  in 
lines  7  and  8.  For  each  pipe  section  present  in  the  network  as  shown  in 
the  order  given  in  lines  7  and  8,  two  lines  of  information  are  required. 

The  first  specifies  the  length  of  each  pipe  in  that  section  in  the  order 

of  direction  of  flow;  the  format  is  10F10.4.  The  second  line  contains 
the  diameter  data  in  the  same  order  and  format  as  the  length  data.  Even 
though  the  format  of  each  line  is  10F10.4,  only  data  for  the  number  of 
lengths  in  that  section  are  required.  If  a  section  is  not  present  no 
lines  of  data  are  input  for  it.  If  a  throttling  device  is  to  be 
simulated,  a  value  of  -99-0  is  input  for  length  and  diameter.  (Area  and 
discharge  coefficient  data  for  the  throttling  device  are  inserted  into 
the  throttling  device  subroutines  as  discussed  in  the  next  section.) 

Last  line  -  Contains  the  values  of  specific  heat  ratio  and  gas 
constant  (ft-lbf )/(lbm-°R)  in  that  order  and  in  format  2F10.4. 


THROTTLING  DEVICE  AND  TANK  CONDITIONS 


Throttling  device  and  tank  condition  data  are  still  required  to 
completely  define  the  network.  Since  the  respective  subroutines  are 
based  on  equations  that  are  peculiar  to  a  particular  type  of  throttling 
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device  or  tank  condition,  the  user  must  edit  the  program  to  change 
equations  and  constants  for  discharge  coefficient,  area,  mass  flow  rate, 
expansion  factor,  etc.  The  subroutine  source  listings  contain  comment 
statements  to  help  the  user. 

Throttling  device  calculations  are  performed  by  subroutine  FVALVE  for 
pipe  sections  1  through  4  and  subroutine  BKVALVE  for  pipe  sections  5 
through  14.  In  these  subroutines  there  is  a  separate  calculation  for 
each  pipe  section  so  that  different  throttling  device  characteristics  can 
be  defined.  The  mass  flow  rate  and  expansion  factor  equations  that  are 
presently  used  in  the  subroutines  are  shown  in  the  "Theory"  section  of 
this  report.  Should  these  not  apply  to  a  particular  throttling  device, 
then  other  equations  must  be  substituted.  In  any  case,  discharge 
coefficient  and  related  area  must  be  inserted  in  the  form  of  a  constant 
or  a  function  of  time  (time  is  transmitted  to  the  subroutine  through  its 
argument) . 

Subroutine  CTANK  calculates  tank  conditions  at  the  end  of  each  time 
step  using  the  Runge-Kutta  integration  scheme.  When  a  variable  volume 
tank  condition  is  specified  (via  an  integer  value  of  3  on  line  6),  the 
main  program  must  calculate  the  rate  of  change  of  volume  with  respect  to 
time  via  a  user-supplied  differential  equation  prior  to  calling  the  CTANK 
subroutine.  As  explained  earlier,  the  equation  to  define  the  slope  for 
the  variable  volume  must  be  inserted  into  the  main  program.  Similar  to 
the  throttling  device  subroutines,  these  calculations  use  separate 
calculation  sections  in  the  main  program  for  each  tank.  The  CTANK 
subroutine  assumes  an  isothermal  thermodynamic  process  between  time 
3teps,  but  can  be  changed  by  the  user  to  include  other  processes  as 
required . 

This  concludes  all  requirements  for  defining  the  network. 


PROGRAM  EXECUTION 


The  compressible  flow  computer  program  is  presently  configured  to  run 
on  a  Digital  Equipment  Corporation  VAX  11/780  computer  with  the  user 
providing  setup  data  via  interactive  terminal.  Prior  to  any 
calculations,  the  program  asks  the  user  for  the  input  and  output  data 
file  names,  the  time  increment  interval  in  number  of  time  steps  at  which 

results  are  to  be  written  to  the  output  file,  the  integration  time  step 

size  (s),  the  maximum  run  time  (s),  and  the  initial  estimate  of  mass  flow 
rate  (lbm/s)  from  tank  1.  This  is  all  the  information  required  to 
execute  the  program.  Note  that  if  the  user  sets  integration  time  step 
size  equal  to  the  total  run  time,  then  the  program  provides  a 
steady-state  network  analysis. 

Program  termination  occurs  when  one  of  three  criteria  are  met:  when 

the  maximum  specified  run  time  has  been  reached,  when  the  pressure  in  the 

supply  tank(s)  is  less  than  or  equal  to  that  in  the  receiving  tank(s),  or 
when  flow  from  tank  1  has  stopped.  As  mentioned  in  the  "Theory"  section, 
tank  1  is  the  main  supply  tank  and  controls  the  calculations;  hence  it 
must  always  produce  a  flow  in  the  downstream  direction. 


OUTPUT 

Appendix  C  contains  the  output  file  which  resulted  from  executing  the 
program  using  the  sample  input  file  of  appendix  B.  Note  that  all  input 
data  contained  in  the  input  data  file  are  initially  written  to  the  output 
file,  permitting  verification  of  input  data.  Results  of  program 
calculations  at  each  time  step  for  which  results  were  to  be  written  are 
listed  next.  At  each  output  time  step,  the  following  data  are  written  to 
the  output  file: 

1.  Time  averaged  mass  flow  rate  (lbm/s)  through-  each  section 


Time  averaged  stagnation  pressure  at  the  end  of  each  pipe  section 
Time  averaged  Mach  numbers  at  the  end  of  each  pipe  section 
Tajik  pressures  at  the  end  of  the  time  step 
Tank  volumes  at  the  end  of  the  time  step. 


SUMMARY  AND  CONCLUSIONS 

A  method  of  calculating  the  transient  flow  of  compressible  fluids 
through  complex  piping  networks  has  been  developed,  and  a  computer 
program  implementing  this  method  has  been  written.  This  report  has 
documented  the  theory  and  logic  upon  which  the  solution  method  and 
computer  program  are  based.  It  also  serves  as  an  introduction  to  the 
p.ogram,  which  was  designed  to  be  used  by  persons  having  minimum 
knowledge  of  compressible  flow  theory. 

The  program  i3  shown  to  be  sufficiently  general  to  handle  a  wide 
variety  of  complex  piping  systems.  Supply  and  receiving  tanks  can  each 
be  defined  as  having  either  constant  volume,  constant  pressure,  or 
variable  volume;  flow  can  be  either  steady  or  unsteady;  flow-limiting 
conditions  resulting  from  wall  friction,  changes  in  pipe  diameter,  and 
throttling  devices  can  be  simulated;  and  the  number,  length,  diameter, 
and  configuration  of  pipe  sections  can  be  easily  varied.  Numerical 
techniques  such  as  Runge-Kutta  integration  and  Newton-Raphson  iteration 
have  been  incorporated  in  the  program  wherever  possible.  The  program, 
besides  being  versatile  and  adaptable  to  a  wide  variety  of  piping 
networks,  provides  rapid  solution  convergence. 
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APPENDIX  A 
FT.OWCHAHTS 


wain  program 


CALL  SUBROUTIN*  RATIQ  TO 
A04UST  FL0W<1) 


MILS 

TEST 


DIMENSION  all  VARIASLES 

u  “  1 

RCAO  INPUT  DATA  PILE  ] 

DETERMINE  INITIAL  MASS  IN  EACH  TANK 

U_  | 

ESTIMATE  FLOW!!) 

t 

CALL  SUBROUTINE  CTANK  TO  PINO 
PRESSURES  ANO  PLOW  RATES  PROM 
SUPPLY  TANKS:  CALCULATE  PRESSURE  AT 
COMMON  POINT  'JUNCTION  B> 

I 

▼ 

CALL  SUBROUTINE  REVERSE  TO  PINO 
PRESSURES  ANO  PLOW  RATES  IN 
RECEIVING  TANKS  ANO  ALL  PIPES  BACK  TO 

junction  a 

i 

COMPARE  PRESSURES  At  JUNCTION  S 
QETA1NE0  PROM  SUBROUTINES  CTANK 

ANO  REVERSE.  (THIS  IS  THE  OVERALL 
CONVERGENCE  TEST) 

PASSES  |  TEST 

OVERALL  PASS  IS  COMPLETE;  COPPECT 
FLOWd)  IS  KNOWN;  CORRECT  MASS  PLOW 
RATE  DISTRIBUTION  ANO  JUNCTION 
STAGNATION  PRESSURES  FOUND 

i 

WRIT*  OUT  RtSULTS 

1 

r 

HAS  MAXIMUM  TIM*  SEEN  EXCSEDCO  T 

NO 


AS  A  FIRST  EST1MAT*  US* 
0L0  FLOW(I) 


CAU  SUBROUTINE  RUNOE  TO 
UPOATS  tank  PRESSURES 
ANO  VOLUMES 


res 


STOP 


A-i 


$ 


£ 

“l  f 


SUBROUTINE  CTANK 


«t  .*•  ~ *.•+'*' *■*+*•  ****** *'*,•* 


u*.***^ jfi\A^MJL».ki.mj*>L>i»:iilu>t  >.>,t  «  i.i  i.>r4 1'4  ^♦j  lS  ivk'i 


ADJUST  FLOW  (2)  VIA 
SUBROUTINE  RATIO 


ADJUST  FLOW  (3)  VIA 
SUBROUTINE  RATIO 


wri:  THIS  FLOWCHART  REPRESENTS  THE 
MOST  COMPLEX  CASE  IN  WHICH  FIFE 
SECTIONS  1  THROUGH  3  ARE  ALL 
FRESENT.  CASES  WHERE  NOT  ALL 
SECTIONS  ARE  PRESENT  ARE 
AUTOMATICALLY  HANDLED  BY  THE 
PROGRAM.  WHICH  COMPENSATES 
FOR  THE  MISSING  SECTIONS  ANO 
VERIFIES  FLOW  CONTINUITY 
CHOKING  LOGIC.  WHICH  IS  0ETA1LE0 
IN  THE  BOOT  OF  THIS  REPORT  AHO 
IN  THE  COMPUTER  PROORAM  ITSELF. 
HAS  BEEN  SUPPRESSED  IN  THIS 
FLOWCHART. 


DIMENSION  ANO  INITIALIZE 


CALCULATE  FROM  TANK  1  TO  JUNCTION  A 
VIA  FIFE  SECTION  1 


CALCULATE  FROM  TANK  2  TO  JUNCTION  A 
VIA  PIPE  SECTION  2 


CHECK  STAGNATION  PRESSURES  FROM 
PIPE  SECTIONS  t  AN 0  2  AT  JUNCTION  A 


00  PRESSURES  CONVERGE  ? 


CALCULATE  FROM  JUNCTION  A  TO 
JUNCTION  B  VIA  PIPE  SECTION  4  USING 
FLOW  (4)  -  FLOW  (7)  ♦  FLOW  (2> 


CALCULATE  FROM  TANK  3  TO  JUNCTION  » 
VIA  PIPE  SECTION  3 


CHECK  STAGNATION  PRESSURES  FROM 
PIPE  SECTIONS  3  ANO  4  AT  JUNCTION  B 


00  PRESSURES  CONVERQC  ? 


POCTANK  ■  STAGNATION  PRESSURE  AT 

junction  b 
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SUBROUTINE  3XVALVS 
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vis 
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SUBROUTINE  REVERSE  (CONT‘0  ) 


rlQW(7)  •  FLCWCTANK- 
_ FLQW(fl) 


FLOWl#)  -  RATIO  -  F1.0WC7ANK 


ADJUST  RATIOS 
WITH  SUBROUTINE  ratio 


CALCULATE  FROM  JUNCTION  g  TO 
JUNCTION  0  VIA  PIPt  SECTION  8 


RATIOS  -  FLOWIBVFLOWCTANK 


CALCULATE  FROM  JUNCTION  0  TO 
JUNCTION  C  VIA  PIPE  SECTION  6 


CALCULATE  FROM  JUNCTION  0  TO 
JUNCTION  C  VIA  PIPE  SECTION  7 


CHECK  STAGNATION  PRESSURES  FROM 
PIPE  SECTIONS  S  ANO  7  AT  JUNCTION  C 


DO  PRESSURES  CONVERGE  ? 


I  YES 


CALCULATE  FROM  JUNCTION  C  TO 
JUNCTION  3  VIA  PIPE  SECTION  3 


PORCVER3E  ■  STAGNATION  PRESSURE 
Of  PIPE  SZCVQN  S  AT  JUNCTION  S 


NOTE;  THIS  FLOWCHART  REPRESENTS  THE 
MOST  COMPLEX  CASE  IN  WHICH 
PIPE  SECTIONS  S  THROUGH  t4  ARC 
ALL  PRESENT.  CASES  WHERE  NOT 
ALL  PIPE  SECTIONS  ARC  PRESENT  ARC 
AUTOMATICALLY  HANOLEO  SY  THE 
PROGRAM.  WHICH  COMPENSATES 
FOR  THE  MISSING  SECTIONS  ANO 
VERIFIES  FLOW  CONTINUITY. 
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SUBROUTINE  PRES 


sufer-choked  inlet 


UNCHOKEO  INLET: 

FLO  -  MOOT 

CALCULATE  MACH1  BASED 
ON  BA  no 


t  r  .  >  i  ,  .  .  »-  l.'*  . »  . 


WW.AAAf.'itWiC’iVlVl 


CALCULATE  RATIO  .  A/A* 
BASED  ON  MOOT 


RATIO.  LE.  1.3  1 

~~  |  TES 

PLOWMAX  «  MOOT 


CALCULATE  RATIO  •  A/A* 
BASED  ON  PLOWMAX 


RATIO.  LT.  1.0  ? 

[  YES 

FLOWMAX  .  0-395*  FLOW  MAX 


|  MACH1  »  1.0  j 

_ i _ 

EVEN  IF  RATIO.  £0.1.0  AT  THE  INLET 
FRICTION  WILL  CAUSE  SUPER-CHOKED 
CONDITION:  SWITCH  »  3;  FLO  •  PLOWMAX 


ESTIMATE  FRICTION  FACTOR 
FROM  MACH1  AND  FLO 


CALCULATE  AN  INLET  MACH  NUMBER 
(MACH1MAXI  THAT  WILL  JUST  CAUSE  THE 
FIFE  TO  CHOKE  AT  THE  EXIT 


USE  MACH1MAX  TO  CALCULATE  A 
CORRESPONDING  FRICTION  FACTOR 
(FI  MAX) 


MAS  F1MAX  CONVERGED  ? 


1 "  AiWlW 


SUBROUTINE  PRES  iCQNTQ  ) 


YES 


I*  (MOOT.  GT  FLCWMAX)  SWITCH  »  : 


.1  i.«. 
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APPENDIX  3 

SAMPLE  INPUT  DATA  FIL 


3AM»Lr  N5T*CRK 
540.0 

2033.0,2083.0,2033.0,100.0,140.0,100.0,140 

320. . 320. .320. .200. .200.  ,200. ,200. 

1,1, 1,1, 1,1,1 

2, 2, 2, 2, 2, 2, 2 

3. 3. 3. 2. 2. 2. 2 

2. 2. 2. 2. 2. 2. 2 

12.  ,13.  ,14. 

.Ill, .  08  0  , .0  70 

1 2  .  ,  1 4  . ,  1 2  . 

.  0  35  , .  081  , .  077 

14.  ,  16. , 15. 

.  0  7  5  ,.  03  1  ,  .0  77 

13  .  ,-99 

.  0  39  ,  .  111 

15.  ,  16. 

.125, .105 

14 . .  1 6. 

.  0  3  9  ,  .  0  8  5 

13.  ,15. 

.155, .166 
15.  ,12. 

.175, .185 

10 . .  13. . 

.035, .038 

12. . 3. 

.079  ,  .083 

13. . -93. 

.0  71  ,  .08  3 
11.  ,11. 

.  0  99  , .  08  3 
12.  ,13. 

.  Ill ,  . 125 

12. .  9. 

.102, .101 
1.4, 53.3 
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APPENDIX  C 

SAMPLE  OUTPUT  DATA  FILE 


$i«ELS  H=T.C3« 


rv^or  ;tr i  t.*sr:« 

i\3UT  5*7*  E[UEN*«S  •  input 

S746E47I0E  TJ»*t'*TU»t  •  5*0.3300  SEJREES  R*H*tNS 

sojctPtc  "t*T  j*t:o  »  i.»ooo 

;i  5  CCiSTInT  •  £3. 5000  Er/3£G«£ES  «»K*I*£ 


'IN*  0‘T» 


TAN* 

• 

l 

STAGNATION 

RRSSSURE 

• 

2750.00 

LB /IN**  2 

VQLSJEE  • 

32.00 

PT**3 

f  a  n* 

« 

z 

S  TAQ. SATIS** 

PRESSURE 

« 

2750.00 

LS/IN**2 

VOLUES  * 

32.00 

*T**3 

**S( 

• 

3 

STAGNATION 

E«  £S  SURE 

■ 

2750.00 

UE/rn»«2 

VSLUEE  * 

32.00 

■  T**3 

r  ana 

• 

li 

stagnation 

PRESSURE 

* 

225.00 

U5/IM**2 

VCLUEE  * 

20.00 

p  r**3 

:  a  ha 

• 

12 

stagnation 

“RES  SURE 

S 

1*0.03 

Ue/IX*»2 

VOLUME  • 

23.00 

e  T**2 

'AS* 

• 

13 

STA  SNATIC-N 

»t*SSU3fi 

s 

22  5. CO 

U3/INPP2 

VCUIES  • 

20.00 

PT**3 

tin  a 

• 

1* 

STAGNATION 

PRESSURE 

* 

1*0.00 

U3/JM»»2 

VOLUEE  • 

23.00 

PT**3 

•45T‘«aa<  •  :•:**£  g 

•  1*5  shct::n  * 

ec-stbt 

i 

•  IP!  a 

1 

lsngth  « 

12.0003- 

*T 

•  :»? 

S5CTI3N 

i 

•  !•?  • 

2 

LENGTH  • 

13.0000 

•T 

<»i  p* 

SSCTIGK 

i 

?r?f  i 

3 

IINGTH  » 

i*.oooo 

•T 

»;9« 

SsCTIGN 

i 

•  ip?  • 

l 

LENGTH  • 

12.0000 

3T 

» r  9  c 

sectign 

z 

ip!  A 

9 

LENGTH  « 

14.0000 

•T 

»:»c 

ssctign 

z 

,»r»5  * 

3 

l  cNGTM  a. 

12.0000 

»T 

•IP? 

SiCTIIN 

3 

UPS  i 

1 

L  ?NGT  N  » 

14.0000 

*T 

SfCTIGN 

3 

•  IPS  1 

2 

LSNGTH  • 

14.0000 

•  T 

a  *og 

SECTION 

3 

PIP?  A 

3 

LENGTH  a 

IS. 0000 

PT 

>Pf 

sjctiin 

4 

PIP?  « 

1 

lsngth  • 

13.0000 

PT 

?t>! 

SJCTI JN 

4 

•  p;  a 

2 

LSNGTH  * 

19.0000 

•T 

ape 

SSCTIIN 

5 

•  IP?  « 

I 

LSNGTH  a 

15.0000 

*T 

apt 

SSCTIGN 

S 

•  IP?  • 

a 

IS  A  CONTROL  06 V IC? 

*:•? 

sscr::* 

5 

•IPS  * 

3 

length  * 

15.0000 

PT 

SSCTIGN 

6 

•  IP?  • 

l 

LENGTH  * 

14.3000 

«T 

*:»£ 

SSCTIIN 

S 

•  IP?  • 

2 

LENGTH  • 

16*0900 

«T 

jp* 

S2CTHN 

7 

•  IP?  i 

1 

LSNGTH  • 

13.0000 

PT 

»!•? 

SSCTI3N 

7 

.*!•?  • 

2 

LSNGTH  a 

1S.0CQ0 

•T 

SSCTIIN 

3 

•IP?  • 

l 

LSNGTH  a 

15.0000 

PT 

»p« 

SSCTIIN 

3 

•  !•£  a 

2 

LSNGTH  a 

12.0300 

PT 

5i»f 

SSCTIGN 

9 

•  IP?  1 

1 

LSNGTH  a 

10. 9000 

•T 

ip* 

SSCTIGN 

9 

•  r»?  « 

2 

LENGTH  a 

1  3.3000 

«T 

ape 

ssctign 

13 

■IP?  • 

I 

length  s 

12.0000 

PT 

ape 

SSCTiin 

13 

•  :•?  • 

2 

LSNGTH  * 

9.3000 

aT 

9  7  9t 

sierras 

1  l 

•ip?  * 

1 

LSNGTH  a 

1  3.3000 

eT 

ip* 

ssctiim 

♦  i 

•IP?  « 

9 

LENGTH  * 

13.0000 

•T 

ape 

ssct:;n 

12 

.»!•?  * 

1 

LSNGTH  a 

11.0300 

PT 

a  I  2? 

sect:  in 

12 

•IP?  a 

2 

LENGTH  * 

11.0300 

•T 

9  *  J« 

pc*:  -N 

l  ! 

•  I  •?  < 

♦ 

. 

length  * 

12-0000 

=  T 

9  *  ae 

ScCt:  ;.s 

»  a 

ape  a 

2 

L  E*lGTH  * 

13.0003 

e  7 

ip* 

SSCTIIN 

: fc 

9  i  9  5  1 

i 

LENGTH  * 

12.C3Q0 

*  T 

a  -  aa 

aja«  1 

a 

LENGTH  a 

9  •  3  0  C  3 

x  T 

or*»eTE< 
c:*»€TE3 
3I4EETE5 
Cli-ETE* 
C£1*<ETE* 
::*<etei< 
oi»««te.< 
ounete^ 
oueete* 
CI»«ETE* 
0  li  “E  Tf  5 
0  t*««TS.< 

0I*«ETE.< 

CHESTER 

Ot**ETE< 

:i**ETE3 

::*eete* 

CU^STER 

CUEETER 

o  ireetsr 
Chester 
:t  j“ET53 
CtiEfTER 
0I*e*TE< 

o:i*st«* 

3:j»ST£.R 

C:i"E*E< 

::***ts» 

::3“ET*3 

ill';’!* 


0.1110  «T 
0.0500  *T 
0.0700  -T 
0.0550  «T 
O.OEIO  ‘T 
0.0770  =T 
0.0750  *7 
0.0510  «T 

o.ott:  et 

0.0990  T 

o.mo  ■’ 

0.1250  ET 

0.0590  ET 
0.C993  *T 
0.0550  *7 
0.1553  ET 
0.1440  *7 
0.1753  ET 
0.1350  *T 
0.0550  «T 
0.09«0  *T 
0.0790  *T 
0.0133  *T 
0.0710  JT 
0.333  3  * T 
0.3453  “7 
0.0740  *T 
3.0-53  =T 
0.1253  e 7 
3.1320  *7 
2.1310  c  r 


C-l 


si 
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SAMPLE  OUTPUT  DATA  FILE 


F 


V 


$ 


c 

i 


r 


V 

•> 

s 

\ 

\ 


r 

l 

8 

s 

V 


$A«»L*  NE7.C3A 


:n»ut  :m  listing 

£.\»ur  oat*  ailena*'  *  in’ut 

stagnation  te.aeaatu.e  »  5*0.0300  oegaees  oankins 

s^eciaic  »eat  9*r:o  »  i.»ooo 

GAS  CC  tSTi-iT  »  50.  3000  5T/3EG9EES  «»N«IN« 


tan*  o*T» 


TANK 

t 

l 

stagnation 

99SSSL9E 

a 

2750.00 

t,E/IN»*2 

VOLUME  * 

32.00 

FT**  3 

r  ANA 

• 

2 

stagnation 

o«  55  SUAE 

a 

2750.00 

LS/IN..2 

V3Lj««  ■ 

32.00 

6Tm  3 

’ANA 

• 

3 

STAGNATION 

’«  £S  5U  9  g 

« 

2750.00 

LS/!n»»2 

V2LJFE  * 

32.30 

«  T*»3 

r  INK 

• 

11 

stagnation 

99  ESSU«E 

• 

235.00 

LB/IN*»2 

VCUUM6  * 

20.00 

F  T**3 

’INK 

« 

l; 

S'agnat;;n 

®6  SS  SL’E 

a 

1*0.31 

UE/In«*2 

vgwl'ns  - 

23.39 

*fi«  3 

’ANA 

t 

13 

stagnation 

*9  SSSL5E 

a 

225.00 

LS/IN*«2 

VOLUM?  « 

23.00 

FT *■*  3 

tank 

i 

14 

STAGNATION 

99JSSU9E 

a 

1*0.00 

L  5 / 1 N**2 

V3LJHE  ■ 

20.00 

CT  *«3 

4= 'aC*  <  Sl.aiNi  a 

®:>E  s£ct:;h  « 

IC-ETOr 

l 

?t?«  • 

I 

length  • 

12.0000 

«T 

3!  if 

SSCTI2H 

1 

’I’E  t 

2 

LENGTH  * 

13.0000 

FT 
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RESULTS  AT  END  05  ®IM  SECTIONS  l  TO  l* 
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